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Abstract 

We present the theory of analytical Hartree-Fock gradients for periodic sys- 
tems as implemented in the code CRYSTAL. We demonstrate how derivatives 
of the integrals can be computed with the McMurchie-Davidson algorithm. 
Highly accurate gradients with respect to nuclear coordinates are obtained 
for systems periodic in 0,1,2 or 3 dimensions. 



Typeset using REVTpjX 



I. INTRODUCTION 



The determination of equilibrium structure is one of the most important targets in elec- 
tronic structure calculations. Analytical gradients provide an important tool to facilitate 
this and therefore the implementation of analytical gradients has become an important part 
of modern codes. Although most solid state calculations are nowadays performed within the 
framework of density functional theory, Hartree-Fock theory can serve as a useful starting 
point for a correlation treatment. In the field of quantum chemistry, a Hartree-Fock so- 
lution is necessary to make a wavefunction based correlation scheme such as, for example, 
the coupled-cluster approach, applicable. Therefore, the determination of a Hartree-Fock 
solution is often an important target. 

The calculation of Hartree-Fock gradients was pioneered by Pulay who performed the 
first implementation for multicenter basis sets0. It should be mentioned, that the theory had 
already been derived earlier independently!. Analytical gradients have become an important 
area in quantum chemistry and several review articles have been published^!. 

Significant work has already been performed for one-dimensional systems: formulas for 
analytic gradients, with respect to nuclear coordinates as well as with respect to the lattice 
vector, have been derived and implemented in a periodic code@; and the theory has been 
extended to metallic systems!! Further progress has been the derivation and implementation 
of formulas for MP2 energyi0 and gradients^, as well as gradients on the density functional 
leveli. Even formulas for second derivatives have meanwhile been coded^l. Recently, a 
scheme for an accurate treatment of long-range Coulombic effects in Hartree-Fock gradients 
has been presente and a new implementation of density functional energy and gradients 
for periodic systems has been demonstated to be highly efficient and accurate^!. 

In this article, we report on an implementation of Hartree-Fock gradients with respect 
to nuclear coordinates in a general periodic code (the CRYSTAL§Hll package), which is to 
the best of our knowledge the first implementation for the case of 2- and 3- dimensional 
periodicity. 
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The article is structured as follows: in section 0, the basis functions and Hartree-Fock 
equations are given. The calculation of integrals which relies on the McMurchie-Davidson 
algorithm and the calculation of gradients of the integrals is explained in sections [TIT] and 
[TV| . The total energy as calculated by the CRYSTAL code is given in section |V[ In section 
[VT| , we explain the calculation of forces and possible sources of error. Finally, in section |V 11 
we illustrate the accuracy of the gradients with some examples. 

II. BASIS FUNCTION AND HARTREE-FOCK EQUATIONS 

In this section, we summarize the basis functions used in the CRYSTAL code and give 
the structure of the Hartree-Fock equations. 



A. Basis functions 

Unnormalized spherical Gaussian type functions (SGTF) in a polar coordinate system 
characterized by the set of variables (|f|, </?), and centered at A, are defined as 

S(a, r — A, n, I, m) = \f — A\ + P) m '(cos $) exp(im<yj) exp(— a\r — A\ ) (1) 

with p) m ' being the associated Legendre function. In the context of the McMurchie- 
Davidson algorithm, Hermite Gaussian type functions are necessary which are defined as: 




Real spherical Gaussian type functions are defined as 

R(a,r — A, n, I, 0) = S(a,f— A, n, I, 0) 
R(a,r — A, n, I, |m|) = Re S(a, r — A, n, I, \m\) 
R(a,r — A,n,l,—\m\) = Im S(a,r — A,n,l, \m\) 

CRYSTAL uses real spherical Gaussian type functions, which are in the following denoted 
as 4>,j,{f — Afj) = N^R(oi, A/t, n, I, m), with the normalization N^. \x is an index enumerating 
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the basis functions in the reference cell (e.g. the primitive unit cell). Although the code 
allows only the use of SGTFs with n = 0, in the process of the evaluation of molecular 
integrals, SGTFs with n^O are usecH 

B. Hartree-Fock equations 

The Hartree-Fock treatment of periodic systems^ is briefly repeated in this section. 
We assume that orbitals are doubly occupied and work within the restricted Hartree-Fock 
formalism. The crystalline orbitals are linear combinations of Bloch functions 

%(r,k)^J2 a ^)Mr,k) (3) 

which are expanded in terms of real spherical Gaussian type functions with fixed con- 
traction coefficients df 

</v(f, k) = E djR(a j7 f-A^-^n, I, m)e ik ° (4) 

g J 

The sum over g is over all direct lattice vectors. The Hartree-Fock-Roothaan equations 
have a structure similar to the molecular case and have to be solved on a set of points in 
the reciprocal lattice: 

E F^(k)a vi (k) = S^(k)a ui {k)ei(k) (5) 

V V 

The Fock matrix F is given in detail in section of this article, the overlap matrix S in 



section pj. The spin-free density matrix in reciprocal space is defined as 

P^{k) = 2 E tv(£K(£)6(e F - ti{k)) (6) 

i 

with the Fermi energy ep and the Heaviside function O; i is an index enumerating the 
eigenvalues. The density matrix in real space is obtained by Fourier transformation. 
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III. CALCULATION OF INTEGRALS 



A. Types of integrals 

In this section we summarize the appearing types of integrals. This is done with the 
assumption that the basis functions, are real. 

1. Overlap integral 
The basic integral is the overlap integral: 

3 wiv& = J A* - 9i)Mr~ K - <? 2 )d 3 r (7) 

The integration is over the whole space, i.e. x from — oo to +00 and similarly for y and 
z. Exploiting translational invariance we can rewrite this as 

Sfrg = J Mr- A,)Mr-A u - g)dh (8) 

with g = g 2 - g[. 

Further integrals appearing are: 

2. Kinetic energy integrals 
Tfrs = J M?- A»)(-l&r)Mr- K - ^)d 3 r (9) 



3. Nuclear attraction integrals 
Nfrg = ~ E Za I Mr- A,)A(f- A a )Mr- K - g)c\\ (10) 

a J 

where A is defined as the Coulomb potential in the molecular the Euler- 

MacLaurin potential for systems periodic in one dimensionil'ii, as Parry's potential^ for 
systems periodic in two dimensions, and as the Ewald potential for systems periodic in three 
dimensionsil'ii. The a summation runs over all nuclei of the primitive unit cell. 



4- Electron- electron Coulomb interaction integrals 



A bielectronic integral can be defined as 
B ^ = f <i>y(r- A^ v {r- A v - g)<p T {.r' - A T - n)<j) a (r' - A a - n - jj) 3 3 , 

fAOvgrnan+h J lj? ^»/| 

(11) 

In the context of periodic systems, it is necessary to perform summation over all lattice 
vectors g, h, n. We define a Coulomb integral as follows 

_ /• <j>u(r- A a )<p u (r- A v - g)(j> T {r' - A T - w)<ft g (r ' -A a -n-h) 3 3 , _ 

^uOugrOah ~ 2^ J \f_ f/\ CL r Q r 

pen 

YB* (12) 

n 

The penetration depth pen is defined as those terms for which 

J [<?« n ,f- A,)g{aT\r- A u - g)\ l ' 2 g{oC n , r- A T - n)d 3 r > l(r' TOi2 (13) 
holds, with g(a™ in ,r) = (^L) 3 / 4 exp(-a™ n |^ 2 ). 

g(a™ n , r) means the lowest exponent of all Gaussians centered at A a . For these integrals, 
the Coulomb interaction is evaluated without approximation. All the other integrals are 
evaluated with a multipolar expansion. "JTOL2" is a tolerance which can be chosen by the 
user of the code. This criterion introduces an asymmetry in the energy expression: a given 
bielectronic integral might be evaluated in different ways for B » - - - , ? and B - - . ? * 

00 J nwgrnon+h rnan+h^iUug 

Avoiding this would however be very inefficient: To keep the symmetry, the lattice sum 
would have to be further broken down into pieces which are evaluated exactly and other 
pieces which are approximated — this would require a much higher computational effort 
and more disk storage. The simpler criterion in equation [13] minimizes the effort, and, when 
ITOL2 is chosen sufficiently large, the violation of the symmetry is negligible. We will 
illustrate this with some examples in section [VT1. 
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5. Electron- electron exchange interaction integrals 
v _ /' 4>i*(f- An)<f>r(f- A T - n)(j) v {r' -A v - g)(j) a (f' -A a -n-h) 33f 

X ^vgrOah ~ 2^ J \f-f'\ 1 ~~ 

^ ^ ^ fiOrnvgcrri+h (~^) 
n 

For an individual exchange integral, 

^ fiOrnvgcrri+h ^Tn/iOcrn+hug ^) 

should hold. However, for efficiency reasons, two different thresholds have been 
introduced@ which leads to another possible asymmetry: an exchange integral is discarded 
when the pseudooverlap associated with M (r — A^) and 4> v {r — A v — g) or the pseudooverlap 
associated with T (r — A r — n) and <f) a (r — A a — n — h) is smaller than certain thresh- 
olds, io _7TOL4 and 10~ ITOL5 . It is recommended that the threshold ITOL5 associated with 
4> T (f— A T — n) and (j) a (f — A a — n — h) should be higher than IT OLA. This, however, will 



lead to a violation of equation [15] and therefore another asymmetry in the energy expression. 
A further cutoff parameter, ITOL3, selects the exchange integrals symmetrically: exchange 
integrals are also neglected if the overlap of 0^(r — A^) with T (r — A T — n) or the overlap 
of (j) v (f— A v — g) with <fi a (r — A a — n — h) is lower as 10 _7TOL3 . This is a symmetric cutoff, 
and therefore should not lead to an inaccuracy in the forces. 

6. Multipolar integrals 

The charge distribution is approximated with the help of a multipolar expansion up to 
order L. 

r,r( P ; A c ) = J p(f)Xr(r-- I c )d 3 r (16) 
with X™ being regular solid harmonics^ and the charge p(r) defined as 

p(f) = - £ P u ^Mr- A»)Mr~ A v - g) (17) 
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7. Field integrals 



The electrostatic potential is approximated with an expansion up to the maximum quan- 
tum number L. 



pen 

IfiQi/ge 



/r Pel i -i 

<^(f - AJMr- K - g) A(f- A c ) - E — 
- \r — A 



\f—A r — n\ 



d 3 r (18) 



with Z] n {A ( ) being the spherical gradient operator, Z'\ n {A ( ) = ^ 2l 2 1 y A XJ n (A c ) and a 
normalization N™ (Ref. p6| ). 



8. Spheropole 

This term arises because the charge distribution is approximated by a model charge 
distribution in the long range. However, the use of the Ewald potential instead of the 
Coulomb potential requires a correction in the three-dimensional caseB 

Q = HQc = H^j (PeW - P r dcl (rlW\ 2 dh (19) 

with 

Pc(r) = - E E Pu^oMr- A,)Mr- X - g) (20) 

and 

pr dcl (r) = e E c(p c ; A^rA, #o (2i) 

Z=0 m=-l 

and 

<H£, f) = Urn Z™(A c )Hoc, f - A c , 0, 0, 0) (22) 
c is chosen as the set of basis functions sited at center A r . 
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B. McMurchie-Davidson algorithm 



In this section, we indicate how the integrals are evaluated. This is done with the 
McMurchie-Davidson@ algorithm. In this formalism, the product of two Gaussian type 
functions is expanded at an intermediate center in terms of Hermite Gaussian type functions. 

S(a, r — A,n, I, m)S(j3,r — B, h, I, m) = 

E(n, I, m, h, I, m, t, u, v)A(j,r — P, t, u, v) (23) 

t,u,v 

with 7 = a + 3 and P = aA ^~+® ■ E also depends on a, (3 and the distance B — A; 
however, the dependence on these parameters is suppressed in the notation. 

This makes a very efficient evaluation of integrals feasible. The starting point 
£"(0,0,0,0,0,0,0,0,0) = exp(— ^g\B — A\ 2 ) can be derived from the Gaussian product 
rule§@: 

exp(-a|r - A\ 2 ) exp(-/3|r - B\ 2 ) = exp ^ ~ ^) exp ( ~ ^ a + ^ 

(24) 

The coefficients E can be generated by recursion relations0il. They are zero for the 
case t + u + v>2n + 2n + l + l and for all negative values of t, u or v. 



aA + f)B 
a + 8 



IV. CALCULATION OF DERIVATIVES 

A. Gradients within the McMurchie-Davidson algorithm 

The evaluation of gradients of the integrals is closely related to the evaluation of 
the integrals themselves. All the integrals can be expressed with the help of the E- 
coefficient In the following we show how derivatives of the integrals can be ex- 



pressed in a similar way with Hermite Gaussian type functions. Starting from equation |23|, 
we obtain 



9 



d 



dA x 
d 



(S(a,r — A, n, /, m)S(f3, r — B, n, I, 



m 



OA 



E(n, I, m, n, I, m, t, u, v)A( r y, r — P, t. 



u,v) 



x t,u,v 



d 



( {~QA~ E ( n ' m ' n ' m ' *' M ' V )) A (7,^*- -P, t, u, v) + 
E(n, /, m, n, /, m, t, u, v)A( r y, f — P, t + 1, w, t>) 



t,u,v x 21 

a 



a + (3 

d ex. 
( — — — P(n, m, n, m, t, u, v) + - — — P(n, I, m, h, l,m,t — 1, u, v) ) A (7, f—P,t, u, v)) 



t% KdA x a + P 

G^:{n, /, m, n, Z, m, i, w, u)A(7, r — P, i, w, v) 



(25) 



t,u,v 



Therefore, the gradients can be obtained in a quite similar way as the integrals. Instead 
of the P-coefficients, the coefficients G% (and, for the other derivatives, Gy, Gf, G%, Gy 
and Gf ) have to be used, as defined in equation |25|. We obtain the following relation from 
equation |5]: 



d 



dA : 
a 



G^{n, /, m, h, /, m, t, u, v) = 
■E(n, I, m, n, I, rh, t, u, v) + 



-E(n, I, m, n, l,m,t — 1, u, v) 



(26) 



a + (3 

We could thus derive the G-coefficients from the P-coefficients. However, a more conve- 
nient way would be to have a recursion relation similar to the P-coefficients. Indeed, these 
relations can be obtained in an analogous wayEl. We give the relations here for the case 
of complex spherical Gaussian type functions; a transformation to real spherical Gaussian 
type functions is possible along the lines given in ref. [TR The starting point can be obtained 



from equation |2^ and the definition of P(0, 0, 0, 0, 0, 0, 0, 0, 0) given in section |l IB : 



dA, 



G^(0, 0,0, 0,0, 0,0, 0,0) = 
■P(0, 0,0, 0,0, 0,0, 0,0) + 







a 



a + (3 



£(0,0,0,0,0,0,-1,0,0) 



'a + (3 



(^-^£(0,0,0,0,0,0,0,0,0) 
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and 

g£(o, 0,0, 0,0, 0,1, o,o) = 

d 

E(0, 0,0, 0,0, 0,1, 0,0) + 



dA x 

OL 

E{0, 0,0, 0,0, 0,0, 0,0) 



a + (3 



OL 

——£(0,0,0,0,0,0,0,0,0) (27) 
a + p 

All the other G^(0, 0, 0, 0, 0, 0, t, u, v) are zero. Similarly, we obtain 

G£(0, 0,0, 0,0, 0,0, 1,0) = 
G?(0, 0,0, 0,0, 0,0, 0,1) = 
G£(0, 0,0, 0,0, 0,1, 0,0) 

and 

G^(0, 0,0, 0,0, 0,0, 0,0) = 

a/3 

2 ^5 ( S * - A/W, °' °' °> °' °' °' °' °) 
and 

Gf(0, 0,0, 0,0, 0,0, 0,0) = 

2^-(5 2 - A z )E(0, 0, 0, 0, 0, 0, 0, 0, 0) (28) 
a + p 

Recursion relations for the G-coefficients can be derived using similar arguments as for 
the Z?-coefficients@. There exist recursion relations to generate E{n + 1, /, m, h, I, m, t, u, v), 
E(n,l + l,m,n,l,m,t,u,v), E(n,l + 1,1 + l,n,l,m,t,u,v) and E(n,—l — I,— I — 
1, h, I, m, t, u, v). Recursions are now necessary for G x , Gy and Gf . 

1. Recursion in I and m 

With S(a,r- A, n, I + 1, 1 + 1) = {21 + l)((a; - A x ) + i(y - A y ))S(a,f- A,n,l,l), we 
obtain: 

S(a,f— A, n, I + 1, 1 + 1)S(P, f— B, h, I, rh) = 
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u,v) 



t,u,v 



E(n, I + 1, 1 + 1, n, /, m, t, w, f )A(7, r — P, t 



,u,v) 



t,u,v 



In the case of gradients, we obtain: 

d _ - - - 

-S(a,f— A, n, I + 1, 1 + l)S(P,r — B,n,l,rh) 



dA. 







{21 + 1)— -[ J2 E(n, I, I, n, /, m, t, u, v)((x - A x ) + i(y - A y ))A( T , f - P, t, u, 

OA x \ t ,u,V 

(21 + 1) I —E(n, I, I, h, I, m, t, u, v)A(j, f — P, t, u, v) + 

\ t,u,v 

d ( ~ - 

((x - A x ) + i(y - Ay))-^-yE(n,l,l,n,l,m,t,u,v)A(-f,f- P,t,u,v) 
(21 + 1) I —E(n, /, I, h, I, m, t, u, v)A{ n f, f— P, t, u, v) + 

\ t,u,v 

((x - A x ) + \(y - A y ))G x (n, /, I, n, I, m, t, u, v)A( r y, f- P, t, u, v)j = 
G x (n, I + 1, 1 + 1, n, I, m, t, u, v)A( r y, f — P, t, u, v) 

t,u,v 

By substituting the recursion relation for Hermite polynomials 



(x - P x )A(^,f- P, t, u, v) = — A(7, r - P,t + 1, it, v) + tA(y,r- P,t-l,u,v 

27 



in the penultimate expression of equation BUI we obtain: 



(21 + 1) ^ —E(n, I, I, n, I, m, t, u, v)A(j,r — P, t, u, v) + 



t,u,v 



—A( 1 ,f-P,t+l,u,v) + 

27 

tA(7,f- P,t- l,u,v) + 
(P x - A x )A(<y,r-P,t,u,v) + 
i(—A( 1 ,r-P,t,u + l,v) + 
mA(7, r — P,t,u — 1, f ) + 
(P y - Ay)A(ry,r- P,t,u,v)) )G^(n,l,l,hJ,m,t,u,v) 



From this, we deduce the following relation: 
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G x (n, I + 1, 1 + 1, n, I, m, t, u, v) = 
(21 + 1) ( — E(n, I, I, h, I, m, t, u, v) + 
-^-G x (n, I, I, n, l,rh,t- 1, u, v) + 
(t + 1)G X (n, /, /, n, l,m,t + 1, w, t>) + 
(Pa - A x )Gx (n, I, I, n, I, m, t, u, v) + 
if— G^(n, I, I, h, I, rh,t,u — 1, v) + 
(u + 1)G^(7T,, I, I, n, I, m, t, u + 1, u) + 
(P y — A y )G x (n, I, I, h, I, m, t, u, v)) ) (33) 



In an analogous way, we obtain the following recursion relation: 

G x (n, 1 + 1,-1 — 1, n, /, m, t, u, v) = 

(21 + 1) y — E(n, I, —I, n, I, m, t, u, v) + 
1 

27 



-G x (n, I, —I, n, l,rh,t — 1, u, v) + 



(t + l)G x (n, I, —I, n, l,rh,t+ 1, u, v) + 
(P x - A x )G x (n, I, -I, n, I, rh, t, u, v) - 



i(—G x (n, I, —I, n, /, m,t,u — 1, v ) + 



■JJ_ C A, 

-27 

(u + l)G x (n, I, —I, h, I, m,t,u + 1, v) + 



(P y - A y )G x (n, I, -I, n, I, m, t, u,v))) (34) 



2. Recursion in n 



With S(a,f— A, n + 1, 1, m) — \r — A\ 2 S(a,f — A, n, I, m) the following relation can be 
derived^: 



47 



G x (n + 1, /, m, n, I, m, t, u, v) = 
— (G x (n, I, m, h, l,fh,t — 2, u, v) + 
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G x (n, I, m, h, I, fh,t,u — 2, v) + 

G^{n, I, m, h, I, fh, t,u,v — 2)) + 

-((P x - A x )G^(n,l,m, h,l,fh,t - l,u,v) - 
7 v 

E{n, I, m, h, l,fh,t — 1, u, v) + 
(P y — A y )G^{n, I, m, h, I, fh,t,u — 1, v) + 
(P z — A z )G^{n, I, m, h, I, fh, t,u,v — 1)) + 
(\P - A\ 2 + ( t + U + V + ~ 2) )G^(n, I, m, h, I, fh, t, u, v) - 
2{P X — A x )E{n, I, m, h, I, fh, t, u, v) + 
2((P x -A x ){t+l)G*{n,l,m,h,l,fh,t+l,u,v) - 
(t + l)E(n, I, m, h, l,fh,t+ 1, u, v) + 
{P y — A y ){u + l)G^{n, I, m, h, l,fh, t, u + 1, v) + 
(P z - A z ) (v + l)G^(n, I, m, h, I, fh, t,u,v+ 1)) + 
(t + 2)(t+ l)G^{n,l,m,h,l,fh,t + 2,u,v) + 
[u + 2){u + l)G^{n, I, m, h, I, fh,t,u + 2, v) + 

+ 2){v + l)G^{n,l,m,h,l,fh,t,u,v + 2) (35) 



3. Recursion in I 



Using that S(a,f- A,n,l + l,m) = {{21 + 1)0 - A z )S{a,r- A,n,l,m) - \ f- A\ 2 {1 + 
\m\)S{a,r — A,n,l — l,m))/{l — \m\ + 1) the following relation can be derived^: 



G x {n, I + 1, m, n, I, m, t, u, v) = 
+ y ^ G A {n, I, m, h, I, fh, t,u,v - 1) + 

{P z — A z )G^{n, I, m, h, I, fh, t, u, v) + 
(v + l)G^{n, I, m, h, I, fh, t,u,v + — 
{I + \m\) ((G x (n, I — 1, m, h, I, fh, t — 2, u, v) + 
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G x (n, I — 1, m, h, /, m,t,u — 2, v) + 
G$(n, l-l,m, n, I, m, t,u,v - 2)) /(27) 2 + 
((P x . - AjG^n, / - 1, m, n, l,m,t- 1, u, if) - 
P(n, i — 1, m, fi, I, m, t — 1, u, v) + 
(P y — A y )G x (n, l — l,m, h, I, rh,t,u— 1, v) + 
(P z — A z )G x (n, I — 1, m, h, I, fh, t, u,v — ln/7 + 
(\P - A\ 2 + ( ' t + U + V + ^ )G^(n, l-l,m, h, I, m, t, u, v) - 
2(P X — A x )E(n, l — l,rn, h, I, m, t, u, v) + 
2 ((P x - A x ) (t + l)G*(n, I -l,m, h, l,m,t + l, u, v) - 
(t + l)E{n, I — 1, m, n, l,m,t + 1, u, v) + 
[P y — A y ){u + 1)G X (n, I — 1, m, h, I, rn,t,u + 1, v) + 
(P z - A z )(v + l)G^(n,l- l,m,n,l,rh,t,u,v + 1)) + 
(t + 2)(t+ l)G x (n, l-l,m, n, I, fh,t + 2, u, v) + 
[u + 2){u + l)G x (n, I — 1, m, n, I, m,t,u + 2, v) + 
(v + 2){v + l)G^(n,l - l,m,n,l,rh,t,u,v + 2)^ /(I - \m\ + 1) (36) 

The recursion relations for the G x -coefficients are similar to those for the P-coefficients0, 
with the extra terms arising from the derivative ^-(x — A x ) for the recursion in I and m, 
and from the derivative -^-\r — A\ 2 for the recursion in / and the recursion in n. 

This can similarly be done for Gy and G^. Finally, all G-coefficients are zero for the 
case t + u + v>2n + 2h + 1 + 1 + 1 and for all negative values of t, u or v. 

B. Gradients with respect to other centers 

To obtain the derivatives with respect to center B, the following relation is usedS (note 
R = A-B): 

d d dP x d dR x 



dA r dP r dA r dFLr dA 7 
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ad d 
+ 



7 dP x dR x 



and 



d d dP x d dR x 



dB x 8P x dB x 8R x dB x 
(3d d 



1 dP x dR x 



Therefore, 



dA r dB r dP r 



This means that 



jr^- E ( n ' l ' m ' ^ ™> *> M ' t ') A (7, r- P, t, u, v) 

ati X t,u,v 



d d 



8P X 8A X/ t ^ v 



Applying equation EE we obtain 



G x (n, I, m, n, I, to, t, u, v)A( r f, r — P, t, 



u,v) 



t.u.v 



(37) 



(38) 



Odd 

+ l^ = 7^r ( 39 ) 



J5(n, /, to, n, /, m, t, u, v)A(j, r — P, t, u, v) (40) 



i±7(n, /, to, n, I, to, t, u, v)A( r y, r — P,t + 1, w, i>) 
— ^2 G x (n, /, m, h, I, to, t, u, v)A( r y, r — P, t, u, v) (41) 

t,u,v 

Finally, we conclude: 

(n, Z, m, h, I, m, £, it, v) = E(n, I, m, h, l,m,t — 1, u, v) — G^(n, I, to, n, I, m, £, u, u) (42) 

and similar for derivatives with respect to y and z direction. All the integrals can be 
expressed with the help of the ^-coefficients. Taking the derivative therefore reduces to 
replacing E with the corresponding G^, Gy, etc coefficients. 
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As an example, we show how overlap integral and gradient are obtained. First, the 
overlap is computed as follows: 

/ ^2 E(n, I, m, n, /, m, t, u, v)A( r y, f — P, t, u, v)d 3 r = 



E(n, I, m, n, I, m, 0, 0, 0) (43) 
We have used that 

|A( 7 ,r7,«,t;)d 3 r=^j <WWU (44) 

because of the orthogonality of the Hermite Gaussian type functions (S t o is the Kronecker 
delta). 

The gradient is computed similarly: 
d 



J S(a,f— A, n, I, m)S(/3,r — B, n, I, m)d 3 r = 
/ ^2 Gxi 71 -, h m ; % h ™i t, u -> 1O-M7, r — P, t, u, v)d 



Gl >, i , m , S ,U,0,0,0)g) 3 (45) 

In our implementation, we therefore compute the gradient of the two Gaussians which 
are associated with the integrals, by replacing ^-coefficients with G-coefficients. As a con- 
sequence, if an operator, which might appear in the integral, has a nonvanishing derivative 
(such as, for example, the nuclear attraction), this must be taken into account additionally. 
This derivative with respect to the third center can be obtained by applying translational 
invariance with respect to a simultaneous uniform translation of the three centers. In the 
case of bielectronic integrals, products with two ^-coefficients appear. Obviously, when dif- 
ferentiating, the corresponding rules of differentiating a product must be applied and two 
derivative terms appear, each of them consisting of a product of one set of E— and one set 
of G— coefficients. Finally, the nuclear-nuclear term must be differentiated which is trivial. 
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It is interesting to compare this implementation with that of the Namur groupli-a where 
also gradients within the McMurchie-Davidson algorithm are computed. Whereas our 
scheme computes the derivatives of the two Gaussians appearing in the integral and a possi- 
bly necessary derivative of an operator is obtained by applying translational invariance, the 
alternative implementation^ relies on explicitly computing derivatives of E'-coefficients and 
of the auxiliary function^ Rt,u,v appearing in the integrals. 

V. TOTAL ENERGY 

The total energy consists of kinetic energy, Coulomb energy (nuclear-nuclear repulsion, 
nuclear-electron attraction and electron-electron repulsion), and exchange energy. We as- 
sume that all the orbitals are either empty or doubly occupied. 

A. Kinetic energy 



The kinetic energy of the electrons is obtained as: 

^kinetic = J2 P v ^ I ^(f- A,) ( - lAAMr- K - £)d 3 r 



^vg/iO^fiOug (46) 



B. Exchange energy 

The exchange energy is obtained as: 



B 



exch— el 



irp,, r p - 

^ / t vg^O / j ah+nrn 
9,W h,H,r,a 



I 



<fr^r- A M )0 T (f - A T - n)(j) u {f' -A v - g)<f> a {r' -A a -h-n) 



d 3 r dV 



h,r,a 
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where we have exploited translational invariance of the density matrix with respect to 
direct lattice vectors ft: P ? , -. - = P r 

We can define a Fock operator for the exchange energy which is 

F exch-ei = _\y p x ( 48 ) 

/l,T,(T 



C. Coulomb energy 

Both kinetic energy and exchange energy must converge independently. However, a 
separation of the contributions to the Coulomb energy is not possible: for example, in a one 
dimensional periodic system with lattice constant a, and n being an index enumerating the 
cells, the electron-electron interaction per unit cell would have contributions like: 

OO 1 

(49) 

This sum is divergent (similarly in two and three dimensions). Therefore, in CRYSTAL 
a scheme based on the Ewald method is used to sum the interactionsSS. We only quote 
the results for the individual contributions: 

1. Nuclear-nuclear repulsion: 

e nn = ±J2Z a Z b A(A b -A a ) (50) 

a, b 

2. Nuclear- electron attraction: 

The energy |-Ene ; which is the Ewald energy of the nuclei in the primitive unit cell 
with the all the electrons of all cells, is the same as the energy |-Een, which is the Ewald 
energy of the electrons of the primitive unit cell with all the nuclei in all cells, as long as 
no approximations are introduced^. CRYSTAL uses the following expression as the sum of 
these interactions: 
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7-icoul— nuc \ " r> zpcoul— nuc 

with the Fock matrix p c ° ul - nnc containing the nuclear-electron contributions defined as 
^ nuc = " E Z a J Mr- A,)Mr~ A u - g)A(f- A a )dh (52) 



3. Electron- electron repulsion: 
with the Fock matrix _p c ° ul _7 cl containing the electron-electron contributions defined as 

F Xr l = -Q S fr3 + E PaKroC^grOaH "EE E C(Pc5 A=)M^ c (54) 

lr» c Nra=-1 



D. Total energy 

Finally, the total energy can be expressed as 

^total ^kinetic _|_ ^NN _|_ ^coul— nuc _|_ ^icoul— cl _|_ jjjexch— el 



- F NN _l V P -T - 

g,IJ.,v 

~ E ^oE 2 « / ^(^- 4*)^" £ - 0)^" 4>)d 3 r 

g,fj.,f a 

1 / £ ' 

+ 2 ^ ^W/*o( ~~ Q^^Oug + E PfrhrO^ fiOugrdah ~ E E E ^ (^ c > ^ajM^^ 
g,H,u ' £ T(T c (=0m=-I 

— 7 E ^"W/iO E PahrO-^ fiOugrOcrh (^) 

The Fock operator used has the structure: 

r-itotal /T-i I r-icoul— nuc _i_ r-icoul— cl _i_ riexch— el f^fi^ 

r liOvg ~ 1 fiOvg r tiOug ^ %®vg ^ %ttvg ^> 

We note that this expression for the Fock operator would be exact if we could guar- 
antee that the penetration depth and screening was symmetric. This would require that 
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^aOv§TOah = ^ rQah^ovg should always hold. This, however, as aforementioned, cannot be 
guaranteed because the truncation is not necessarily symmetric. In addition, the screening 
of the exchange interaction is not necessarily symmetric. Therefore, an inaccuracy in the 
Fock operator will show up which will be stronger the more asymmetric the truncation in 
the energy expression is. 

The total energy can be expressed as 

j-itotal 771NN 1 \ " r> ( T 1 1 77C011I— nuc 1 / Tjicoul— cl , Tjiexch— el\\ ln r 7\ 

E =E + 2- P »gA T rfvg + F ^6ug + 2^ug + %6ug ■>) ^ 



and the Hartree-Fock equations become as in equation |5| In ref. 2f| it was pointed out 
that the quantity QS ^ can be removed from the Fock operator which has been done in 
CRYSTAL. This leads to eigenvalues shifted by Q as we now use the modified equation 

£(F;° tal (£) + QS, v (k))a m (k) =Y,S,*,(k)a vi (k)(e i (k) + Q) (58) 



VI. GRADIENT OF THE TOTAL ENERGY 

The force on the nuclei can be calculated similarly to the molecular caseiS. The deriva- 
tives of all the integrals are necessary, and the derivative of the density matrix is expressed 
with the help of the energy-weighted density matrix. One important assumption is that 

(59) 

holds. Taking the derivative leads, for example, to terms like the following 





8 A 

ur ^* " 9 \W h,T,a 



9,W h,T,<7 

When equation [5£] holds, we rename the indices in the second addend and obtain: 
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^ E ( a J* ^W^O ) E ^crhrd^ \lSugTOah + E ^9>0 X/ ^ctAtO I -* C 



fiOugrOcrh 



(61) 



XOA i 7 h,T,tr h,T,cr ° Ai 

We derived the equation for the force this way although equation does not always 
hold. Therefore, inaccuracies will appear when equation ISD] is strongly violated. The full 
force is obtained as: 



dA 
-VP.- 



9,1*," 



OA 



OA 







+ E ^uguQ E Z a -, 



ugfiO 



OA. 



(f - A^)(j> v {r- A u - g)A{f- A a )d 3 r 



2tt 



d 



^flOuggy E E ^ahrO 
C h,cr,T£c 

r-A T )(f) a (r-A a -h) 



L I 



+ E E / Mr' - AMf 1 -A- h)Xr(r' - A c )d 3 r'5r(A c ,r) 



1=0 m=-l 

h ° hT " dA, 

L I q 

-EEE E 

c i =Qm= _i£ 



r 2 d 3 r 



OA I 
OX 



i T (f - I r )0 CT (f - A a - /I)Xr(r - A c )d 3 r M™^ 



+- V P - - V P - - 

h,r,cr OA 



fiOvgrOcrh 



uu n0vg 



dAi Jbz 



exp(ikg) E 2a "i(*) a «-(*)( e i(*) + - - Q)d 3 k 



(62) 



The last addend is the energy weighted density matrix; the integral is over the first 
Brillouin zone. It is worthwhile mentioning that the factor P V ^S^ V ^ is equal to the number 
of electrons in the unit cell and therefore its derivative with respect to A vanishes. We note 
three important points: 



Equation |62] is correct for the exact solution of the Hartree-Fock equations. Thus, in 
practice, a well converged solution is necessary to achieve accurate forces. 
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The energy- weighted density matrix is fc-dependent. Therefore, the accuracy of the 
forces will become dependent on the number of fc-points. 



The derivation of equation £y assumes that equation [59] holds. The treatment of the 
Coulomb series with finite penetration depth leads to an asymmetry associated with 
ITOL2 in the CRYSTAL code. In addition, in the treatment of the exchange series 
an asymmetry can be introduced if the screening parameters (ITOLA and ITOL5 in 
the CRYSTAL code) are chosen differently. Therefore, the choice of ITOL2, ITOLA 
and ITOL5 will influence the accuracy of the gradients. 



VII. RESULTS FROM TEST CALCULATIONS 

With a few examples, we want to illustrate the accuracy of the analytically computed 
gradients. 

In table |, all the integrals are evaluated without approximation and the analytical 
derivative agrees to five digits with the numerical derivative. As the numerical derivative is 
only accurate up to five digits, this is certainly satisfying. 

In table [TJ the variation in accuracy when penetration depth and overlap criteria are 
altered, is displayed. As described in the article, lowering ITOL2 to low values leads to 
inaccuracies in the gradients. Lowering only one of the parameters ITOLA or ITOL5 also 
leads to inaccuracies, whereas lowering both to the same value gives an analytical gradient 
which is consistent with the numerical gradient — however, as a value of 1 for ITOLA 
and ITOL5 was chosen, energy and force are completely different from calculations with 
reasonable values for ITOLA and ITOL5. The parameter ITOL1, which selects the one- 
electron and Coulomb integrals according to the overlap, can lead to numerical instabilities 
in the energy calculation and inaccurate gradients when chosen too low, and should therefore 
be reasonably high. The parameter ITOL3 does not influence the accuracy of the gradients: 
although when chosen much too low with a value of e.g. 1, numerical and analytical gradient 
still agree. The default values for the JTOL-parameters for the energy calculation are 6, 6, 
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6, 6 and 12. In all calculations performed so far, these default values did not lead to serious 
errors for the gradients. 

Another example (table [TTI| ) is the CO molecule arranged as a single molecule 
("molecule"), as a molecule periodically repeated in one dimension ("polymer"), in two 
dimensions ("slab") and in three dimensions ("bulk"). The forces agree well, and it is 
demonstrated that using stricter real-space truncation parameters improves the agreement. 
The forces seem to be relatively insensitive to the number of sampling points and changing 
their number changed the error in the forces only slightly. 



Finally, in table [TV] we compare analytical and numerical derivatives for MgO when 
moving the oxygen atom in x-direction which would correspond to a longitudinal phonon. 
Again, agreement is to the order of 10~ 5 ^ (default ITOL parameters were used). 



VIII. CONCLUSION 

We presented the theory of analytic Hartree-Fock gradients for periodic systems. This 
has been implemented in the code CRYSTAL which is to the best of our knowledge the first 
implementation of Hartree-Fock gradients in systems periodic in 2 and 3 dimensions. The 
results are in excellent agreement with numerical derivatives. 

Future directions will be the improvement of the efficiency of the code (implementation 
of symmetry and various technical improvements), derivatives with respect to the lattice 
vector, as well as an extension to metallic systems. 
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TABLES 

TABLE I. CO molecule, all tolerances high enough that all the integrals are done without 
approximation. The carbon atom is placed at (0 A, A, A), the oxygen atom at (0.8 A, 0.5 
A, 0.4 A). To calculate the numerical force in x-direction, the x-coordinate of the oxygen atom is 
changed (column 1). The Hartree-Fock energy is displayed in column 2. The numerical force is 
determined from 2 points at (0.8 ±<5) A and displayed in column 3. The analytical derivative for 
x=0.8 A is displayed in column 4. Basis sets of the size [2s2pl<i] were used. 



x-coordinate of oxygen 

A 

0.799 

0.7999 
0.79999 
0.80000 
0.80001 

0.8001 

0.801 



energy 
E h 

-1.107642201574E+02 
-1.107648645627E+02 
-1.107649287000E+02 
-1.107649358230E+02 
-1.107649429453E+02 
-1.107650070153E+02 
-1.107656446886E+02 



numerical derivative 
(x-component) 

E h /a B 
0.376915 
0.376913 
0.376914 



analytical derivative 
(x-component) 
E h /a B 



0.37691274 
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TABLE II. CO molecule, when varying the ITOL parameters. The first and second parameter 
should be reasonably high, the fourth and fifth should both be high or identical to obtain accurate 
forces. 



ITOL values 


numerical derivative 


analytical derivative 




(x-component) 


(x-component) 




E h /a B 


Eh/aB 


6 6 6 6 12 


0.37691 


0.376913 


20 20 20 20 20 


0.37691 


0.376913 


4 20 20 20 20 


0.37691 


0.376912 


3 20 20 20 20 


0.37747 


0.377389 


2 20 20 20 20 


0.38883 


0.388260 


20 2 20 20 20 


0.37684 


0.376821 


20 1 20 20 20 


-7.40676 


-1.246409 


20 20 1 20 20 


-8.09240 


-8.09240 


20 20 20 1 20 


-0.51905 


-3.504654 


20 20 20 20 1 


0.40892 


0.359847 


20 20 20 1 1 


-33.34823 


-33.348229 
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TABLE III. CO molecule, arranged periodic in dimensions ("molecule"), in 1 dimension 
("polymer"), in 2 dimensions ("slab") and in 3 dimensions ("bulk"). We display the accuracy 
of the gradient as a function of the number of sampling points, and as a function of the ITOL 
parameters. 

dimension number of sampling points ITOL parameters numerical derivative analytical derivative 

(x-component) (x-component) 



6 6 6 6 12 



E h /a B 
0.37691 



E h /a B 
0.376913 



5 
5 



6 6 6 6 12 
6 9 6 12 12 



0.37659 
0.37662 



0.376633 
0.376647 



34 



6 6 6 6 12 



0.37630 



0.376335 



3 
3 
3 
3 



260 
260 
260 



6 6 6 6 12 
6 6 6 6 12 
6 8 6 8 14 
8 8 6 8 14 



0.37571 
0.37565 
0.37570 
0.37573 



0.375742 
0.375679 
0.375720 
0.375721 
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TABLE IV. MgO at a lattice constant of 4.21 A. We compare numerical and analytical deriva- 
tives when moving the oxygen ion parallel to the x-direction. Basis sets of the size [3s2p] were 
used. 



displacement of oxygen 

in % of lattice constant 
+ 1 
+2 
+3 



analytical derivative 
(x-component) 
E h /a B 
-0.005923 
-0.012254 
-0.019385 



numerical derivative 
(x-component) 
E h /a B 
-0.00593 
-0.01226 
-0.01940 
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